
% same idea as pool2pVisualStimSessions_v2_0, but this pools the data from
% diving vessels only; this version, instead of grouping based on a global
% vessel ID that is manually assigned, will pool based on location in the
% window
%
%
dowaveletfilter = true;
objType = '25x_v2';
children = dir;
rootdir = cd;
piaTracking = [];
trajectoriesDDD = [];
trajectoriesD = [];
trajectoryinfo = [];
runningdata = [];
eyeDisplacement = [];
eyeRadius = [];
raw = [];
infotable = [];
keyword = 'blinkingCircle';

forbiddenword = 'moving';
for n = 1:length(children)
    if children(n).isdir && children(n).name(1) == 'S'        %iterate through subfolders
        cd([rootdir,'\',children(n).name])
        cursession = str2double(children(n).name(2));
        temp = dir('*groupedDivingResults.mat');
        if isempty(temp)
            continue;
        end
        if strcmp(cd,'\\research.files.med.harvard.edu\neurobio\GU LAB\TWO-PHOTON MICROSCOPY\visualStim\thy1gcamp\013\S2')
            prestimTime = 3.5;
        else
            prestimTime = 4.5;
        end
        try
            temp = dir('*groupedDivingResults.mat');
            if ~isempty(temp)
                regExpindicator{1} = keyword;
                regExpindicator{2} = forbiddenword;
                [a,b,c,d,e,f,~,h,i,j] = poolStimData_diving_v1_4(regExpindicator,[],cursession,objType,prestimTime); 
            else
                [a,b,c,d,e,f,~,h,i,j] = poolStimData_diving_v1_4(keyword,[],cursession,objType,prestimTime);     
            end
        catch
            disp(['unable to do session #',num2str(cursession)])
            continue;
        end
        if prestimTime == 3.5
            a = [nan(size(a,1),10),a];
            d = [nan(size(d,1),10),d];
            e = [nan(size(e,1),10),e];
            h = [nan(size(h,1),10),h];
            i = [nan(size(i,1),10),i];
        end
        infotable = [infotable;j];
        trajectoriesDDD = [trajectoriesDDD;a];
        trajectoryinfo = [trajectoryinfo;b];
        piaTracking = [piaTracking;c];
        trajectoriesD = [trajectoriesD;d];
        runningdata = [runningdata;e];
        eyeDisplacement = [eyeDisplacement;h];
        eyeRadius = [eyeRadius;i];
        raw = [raw;f];
    end
end
cd(rootdir);
if dowaveletfilter
    for n = 1:size(trajectoriesD,1)
        cur = trajectoriesD(n,:);
        missing = isnan(cur);
        if sum(missing) > 0 && sum(missing)<(.25*length(cur))
            cur = resample(cur,1:length(cur));
        elseif sum(missing)>(.25*length(cur))
            continue;
        end
        trajectoriesD(n,:) = wdenoise_NVC_v2(cur);
    end
end

%% calculate baseline for each vessel on a per-loop basis to reduce noise (in the baseline)
baselineframes = 25:45;      %for 4.5s baseline
[DDD,baselinemean,baselinestd,DDDz] = groupDDDcalculate(trajectoriesD,infotable,baselineframes);

%% cluster the data:
numclusters = 4;
clustertype = 'distance';   %use this to just rely on distance

p = trajectoryinfo(:,3:4);
switch clustertype
    case 'distance'
        distanceThreshold = 250;    %in microns
        [Idx,temp] = rangesearch(p,p,distanceThreshold,'SortIndices',true);
        clusterID = zeros(size(Idx));
        for n = 1:length(clusterID)
            curID = max(clusterID(:))+1;
            for n2 = n:length(clusterID)
                if clusterID(n2) == 0 && isequal(sort(Idx{n}),sort(Idx{n2}))
                    clusterID(n2) = curID;
                end
            end
        end
        uid = unique(clusterID);
        infotable.clusterID = clusterID;
    case 'kmeans'
        clusterID = kmeans(p,numclusters);
        infotable.clusterID = clusterID;
        uid = unique(clusterID);
end



%% plot the cluster positions just to verify this worked:
figure
p = trajectoryinfo(:,3:4);
uid = unique(clusterID);

clrs = jet(length(uid));
for n = 1:length(uid)
    scatter(p(clusterID==uid(n),1),p(clusterID==uid(n),2),20,clrs(n,:))
    hold on
end
set(gca,'Ydir','reverse');axis equal
hold off

%% filter out trajectories with too much fluctuation in diameter before stimulus:
bsrange = 10:45;
sloperange = 10:45;
bs = nanstd(trajectoriesD(:,1:45),[],2);        %check for standard deviation in the pre-stim period
maxslope = nan(size(bs));           %check for changes in diameter
for n = 1:size(DDD,1)
    if sum(isnan(DDD(n,sloperange))) > .45*length(sloperange)
        continue;
    end
    temp = smooth(DDD(n,sloperange),10,'rlowess');
    temp = abs(diff(temp));
    temp = temp(1:(length(temp)-5));
    maxslope(n) = max(temp);
end
bsthresh = .35;
slopethresh = .5;
ind = find(bs<bsthresh & maxslope < slopethresh);
DDDfilt = DDD; DDDfilt(ind,:) = nan;

